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Abstract. We present an algorithm for the approximation of a finite 
horizon optimal control problem for advection-diffusion equations. The 
method is based on the coupling between an adaptive POD representa- 
tion of the solution and a Dynamic Programming approximation scheme 
for the corresponding evolutive Hamilton- Jacobi equation. We discuss 
several features regarding the adaptivity of the method, the role of error 
estimate indicators to choose a time subdivision of the problem and the 
computation of the basis functions. Some test problems are presented 
to illustrate the method. 
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1. Introduction 

The approximation of optimal control problems for evolutionary partial dif- 
ferential equations of parabolic and hyperbolic type is a very challenging 
topic with a strong impact on industrial applications. Although there is a 
large number of papers dealing with several aspects of control problems from 
controllability to optimal control, the literature dealing with the numerical 
approximation of such huge problems is rather limited. It is worth to note 
that when dealing with optimal control problems for parabolic equations we 
can exploit the regularity of the solutions, regularity which is lacking for many 
hyperbolic equations. We also recall that the main tools is still given by the 
Pontryagin maximum principle. This is mainly due to the fact that the dis- 
cretization of partial differential equations already involves a large number of 
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variables so that the resulting hnitc dimensional optimization problem easily 
reaches the limits of what one can really compute. The forward-backward 
system which describes Pontryagin's optimality condition is certainly below 
that limit. However just solving that system one is using necessary conditions 
for optimality so, in principle, there is no guarantee that these are optimal 
controls. By this approach for general nonlinear control problems we can ob- 
tain just open-loop control. One notable exception is the linear quadratic 
regulator problem for which we have a closed-loop solution given by the Ric- 
cati equation. This explains why the most popular example for the control 
of evolutive partial differential equations is the control of the heat equation 
subject to a quadratic cost functional. 

In recent years, new tools have been developed to deal with optimal control 
problems in infinite dimension. In particular, new techniques emerged to re- 
duce the number of dimensions in the description of the dynamical system 
or, more in general, of the solution of the problem that one is trying to opti- 
mize. These methods are generally called reduced- order methods and include 
for example the POD (Proper Orthogonal Decomposition) method and re- 
duced basis approximation (see [H]). The general idea for all this method 
is that, when the solution are sufficiently regular, one can represent them 
via Galerkin expansion so that the number of variables involved in this dis- 
cretization will be strongly reduced. In some particular for the heat 
equation, even 5 basis functions will suffice to have a rather accurate POD 
representation of the solution. Having this in mind, it is reasonable to start 
thinking to a different approach based on Dynamic Programming (DP) and 
Hamilton- Jacobi-Bellman equations (HJB). In this new approach we will first 
develop a reduced basis representation of the solution along a reference tra- 
jectory and then use this basis to set-up a control problem in the new space 
of coordinates. The corresponding Hamilton- Jacobi equation will just need 
3-5 variables to represent the state of the system. Moreover, by this method 
one can obtain optimal control in feedback form looking at the gradient of 
the value function. 

However, the solution of HJB equation it is not an easy task from the numeri- 
cal point of view: the analytical solution of the HJB equation are non regular 
(typically, just Lipschitz continuous). Optimal control problems for ODEs 
were solved by Dynamic Programming, both analytically and numerically 
(see [I] for a general presentation of this theory) . From the numerical point of 
view, this approach has been developed for many classical control problems 
obtaining convergence results and a-priori error estimates (0], [Sj and the 
book [S]). Although this approach suffers from the curse-of-dimensionality 
some algorithms in high-dimension are now available ([3] and [2]) and the 
coupling with POD reppresentation techniques will allow to attack by this 
technique optimal control problems in infinite dimension. 
To set this paper into perspective we must say that a first tentative in this 
direction has been made by Kunisch and co-authors in a series of papers 
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[7J [5] for diffusion dominated equations. In particular, in the paper by Ku- 
nisch, Volkwein and Xie [TU] one can see a feedback control approach based 
on coupling between POD basis approximation and HJB equations for the 
viscous Burgers equation. Our contribution here is twofold. The first novelty 
is that we deal with advection-diffusion equations. The solutions to these 
equations exhibit low regularity properties with respect to non degenerate 
diffusion equations so that a rather large number of POD basis functions will 
be required to obtain a good approximation if we want to compute the POD 
basis just once. Naturally, this increases the number of variable in the HJB 
approach and constitutes a is a real bottle-neck. In order to apply the Dy- 
namic Programming approach to this problem we have developed an adaptive 
technique which allows to recompute the POD basis on different sub-intervals 
in order to have always accurate results without an increase of the number 
of basis functions. The second contribution of this paper is the way the sub- 
intervals are determined. In fact, we do not use a simple uniform subdivision 
but rather decide to recompute the POD basis when an error indicator (de- 
tailed in Section 4) is beyond a given threshold. As we will show in the sequel, 
this procedure seems to be rather efficient and accurate to deal with these 
large scale problems. 



2. The POD approximation method for evolutive PDEs 



We briefly describe some important features of the POD approximation, more 
details as well as precise results can be found in the notes by Volkwein [T4] . 
Let us consider a matrix Y £ M. mxn , with rank d < min{m, n}. We will call 
Uj the j— th column of the matrix Y. We are looking for an orthonormal basis 
{^i}i=i ^ m w ith I < n such that the minimum of the following functional 
is reached: 

2 



E 



i=l 



(2.1) 



V = [Vi,...,V n ] £ 
mxn is diagonal, E = 



The solution of this minimization problem is given in the following theorem 

Theorem 1. Let Y = [yi, ■ ■ ■ ,y n ] £ W nxn be a given matrix with rank 
d < min{?7j,n}. Further, let Y — ^>Y,V T be the Singular Value Decompo- 
sition (SVD) of Y, where W = [ipx, ip m ] £ W 
jjnxn are orthogonal matrices and the matrix £ £ 
diag{o~i, . . . , <7 m }. Then, for any £ £ {1, . . . , d} the solution to (2.1) is given 
by the left singular vectors {V'i}|=i> *- e ; by the first £ columns of "J". 

We will call the vectors {^i}| =1 POD basis of rank I. This idea is really 
usefull, in fact we get a solution solving an equation whose dimension is 
decreased with respect to the initial one. Whenever it's possible to compute 
a POD basis of rank £, we get a problem with much smaller dimension of the 
starting one due to the fact £ is properly chosen very small. 
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Let us consider the following ODEs system 

" y(s) = Ay(s) + f(s,y(s)), sg(0,T] 

2/(0) - yo 

where y Q G E m , A G M mxm and / : [0,T] x M m -> 
locally Lipschitz to ensure uniqueness. 



(2.2) 



is continuous and 



The system (2.2 1 can be also interpreted as a semidiscrete problem, where 



the matrix A represents the discretization in space of an elliptic operator, 
say Laplacian for instance. To compute the POD basis functions, first of all 
we have to construct a time grid < ti < . . . < t n = T and we suppose to 
know the solution of (2.2) at given time tj, j = 1, . . . ,N. We call snapshots 



the solution at those fixed times. For the moment we will not deal with the 
problem of selecting the snapshots sequence which is a difficult problem in 
itself, we refer the interested readers to [9]). As soon as we get the snapshots 
sequence, by Theorem [T] we will be able to compute our POD basis, namely, 

Let us suppose we can write the solution in reduced form as 

e i 



substituting this formula into (2.2) we obtain the reduced dynamics 



E vti'Wj = E + f(sy(8)), s g (o,t] 



3=1 



(2.3) 



L j'=i 



We note that our new problem (2.3) is a problem for the £ < m coefficient 
functions yj(s), j — 1, . . . ,£. Thus, the problem is low dimensional and with 
compact notation we get: 

yt(s) = Aty*(s) + F(s,y e (s)) 
2/(0) = yi 



where 



A G 



F = (^,...,F,) T :[0,T]x 



with {A% = (Aip^ipj), 

( Vi \ 

: : [0, T] -> R e 

V v\ J 



Fi(s,y) = (f \ s >£?/j^j ) for s G [0,T] y = (y l} ...y t ) G 
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finally obtaining the representation of yo in 

vl= ; 

In order to apply the POD method to our optimal control problem, the 
number £ of POD basis functions is crucial. In particular we would like to keep 
£ as low as possible still capturing the behaviour of the original dynamics. The 
problem is to define an indicator of the accuracy of our POD approximation. 
A good choice for this indicator is the following ratio 



E <Ti 

i=l 

E <n 

i=l 



(2.4) 



where the a are the the singular value obtained by the SVD. 

As much £ {£) is close to one as much our approximation will be im- 
proved. This is strictly related to the truncation error due to the projection 
of fjj onto the space generated by the orthonormal basis {tp}i =1 , in fact: 

2 



J(Vi 



1=1 



E 



3. An optimal control problem 

We will present this approach for the finite horizon control problem. Consider 
the controlled system 

y( s ) = f(y(s),u(s),s), s£{t,T] , , 

y(t) = xem. n , [ ' 

we will denote by y : [t, T] — > M. n its the solution, by u the control u : [t, T] — > 

R m , f : R n x R m -> R n , s € (t, T] and by 

U = {u: [0, T] U] 

the set of admissible controls where U C lR m is a compact set. Whenever 
we want to emphasize the depence of the solution from the control u we will 
write y(t;u). Assume that there exists a unique solution trajectory for (3.1) 
provided the controls are measurable (a precise statement can be found in 
[T]). For the finite horizon optimal control problem the cost functional will 
be given by 

min J x<t (u) := / L(y(s, u), u(s), s) e - Xs ds + g(y{T)) (3.2) 

where L : R n x R m R is the running cost and A > is the discount factor. 
The goal is to find a state-feedback control law u(t) — $(t/(t), t), in terms of 
the state equation y(t), where $ is the feedback map. To derive optimality 
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conditions we use the well-known dynamic programming principle due to 
Bellman (see pQ). We first define the value function: 

v(x,t) := inf J x t (u) (3.3) 

ueU 

Proposition 3.1 (DPP). For all x£R n andO < r <t then: 

v(x, t) = min { / L(y(s), u(s), s)e~ Xs ds + v(y, t - r) j . (3.4) 
ueu It 



Due to (3.4 1 we can derive the Hamilton- J acobi- Bellman equations (HJB): 
dv 

- — (y, t) = min {L(y, u, t) + Vv(y, t) ■ f(y, u,t)} . (3.5) 
ot ueu 

This is nonlinear partial differential equation of the first order which is hard 
to solve analitically although a general theory of weak solutions is available 
PP. Rather we can solve it numerically by means of a finite differences or 
semi-Lagrangian schemes (see the book [5] for a comprehensive analysis of ap- 
proximation schemes for Hamilton- Jacobi equations) . For a semi-Lagrangian 
discretization one starts by a discrete version of (HJB) by discretizing the 
underlined control problem and then project the semi-discrete scheme on a 
grid obtaining the fully discrete scheme 



,,n+l 



wm[AtL(xi,nAt,u) + I[v n ](xi + AtF(xi,t n ,u))] 
ueu 



with Xi = iAx, t n = nAt, vf v(Xi, t n ) and /[■] is an interpolation operator 
which is necessary to compute the value of v 11 at the point Xi + At F(xi, t n , u) 
(in general, this point will not be a node of the grid). The interested reader 
will find in [5] a detailed presentation of the scheme and a priori error esti- 
mates for its numerical approximation. 

Note that, we also need to compute the minimum in order to get the 
value v" +1 . Since v n is not a smooth function, we compute the minimum by 
means of a minimization method which does not use derivatives (this can be 
done by the Brent algorithm as in [3])- 

As we already told the HJB allows to compute the optimal feedback 
via the value function, but there are two major difficulties: the solution of 
an HJB equation are in general non-smooth and the approximation in high 
dimension is not feasible. The request to solve an HJB in high dimension 
comes up naturally whenever we want to control evolutive PDEs. Just to 
give an idea, if we build a grid in [0, 1] x [0, 1] with a discrete step A:r = 0.01 
we have 10 4 nodes: to solve an HJB in that dimension is simply impossible. 
Fortunatelly, the POD method allows us to obtain reduced models even for 
complex dynamics. Let us focus on the following abstract problem: 

^r(y( s ), <p)n + a (y( s ), <p) = {B{u{s),ip)y,v v^gf 

ds (3.6) 
y(t) =y Q e H, 
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where B : U — > V' is a linear and continuous operator. We assume that a 
space of admissible controls U a d is given in such a way that for each u € U a d 
and yo E H there exists a unique solution y of (3.6 1. V and i? are two Hilbert 
spaces, with (•, we denote the scalar product in H; a : V x V — > M. : is 
symmetric coercive and bilinear. Then, we introduce the cost functional of 
the finite horizon problem 

J yo ,M ■= J L(y(s),u(s), s)e- Xs ds + g(y(T)), 
where L : V x U x [0, T] — > R. The optimal control problem is 

m A n JyoA u ) ( 3 - 7 ) 



subject to the constraint: y g W/ oc (0,T; V) x U solves (3.6) 

with Wi oc {0,T) = f] T>Q W{0,T), where W(0,T) is the standard Sobolev 
space: 

W(0, T) = {<pe i 2 (0, T; V),tft G ^ 2 (0, T; V")}. 



The model reduction approach for an optimal control problem (3.7) is based 



on the Galerkin approximation of dynamic with some informations on the 



controlled dynamic (snapshots). To compute a POD solution for (3.7) we 
make the following ansatz 

t 

y i (x,s) = ^2w i {s)^ i (x). (3.8) 

i=l 

where {V'}f = i is the POD basis computed as in the previous section. 
We introduce mass and stiffness matrix: 

M = ((m y )) € R exe with mij = (i>j,A)H, 

S = ((Sij)) £ R exe with m i:j = a^,^), 

and the control map b : U — > is defined by: 

u b{u) = (b(u)i) £ R* with b{u) l = (Bu, ip,) H . 

The coefficients of the initial condition y £ (0) £ ~R e are determined by Wi(0) = 
(wo)i = (j/0)V')x, 1 < i < £, and the solution of the reduced dynamic 
problem is denoted by w e (s) € Then, the Galerkin approximation is 
given by 

minJ^>) (3.9) 
with u € Uad and w solves the following equation: 

w e (s) = F(w e (s),u(s),s) s > 0, 

(3.10) 

w e (0) = Wq. 



The cost functional is defined: 

Jti M = l T L(w e (s),u(s),s)e- Xs dt + g(w\T)), 
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with w and y e linked to (3.8) and the nonlinear map F : R e x U — > R e is 
given by 

F(w e , u, s) = M^i-Sw^s) + b(u(s))). 
The value function v i , defined for the initial state wo G R £ , 

v e (w e ,t)= inf J l <t {u) 



and w solves (3.9) with the control u and initial condition wq- 



We give an idea how we have computed the intervals for reduced HJB. HJBs 
are defined in R n , but we have restricted our numerically domain which 
is a bounded subset of M". This is justified since y + AtF(y, u) 6 T h for each 
y £ Th and u € U a d- We can chose T/j = [a\, b\\ x [02, ^2] X ■ ■ • M with 
ffli > 02 > . . . > ffl^. How should we compute these intervals [c^, bi]l 
Ideally the intervals should be chosen so that the dynamics contains all the 
components of the controlled trajectory. Moreover, they should be encapsu- 
lated because we expect that their importance should decrease monotonically 
with their index and that our interval lengths decrease quickly. 
Let us suppose to discretize the space control U — {iti, . . . , um} where U is 

symmetric, to be more precise if u € U —u G U. 

1 

Hence, if y e (s) — (y(s), ipil'tpi — w i{ s )i } i^ as a consequence, the coef- 

ficients Wi(s) £ [a%,bi\. We consider the trajectories solution y(s,Uj) such 
that the control is constant u(s) = uj for each tj, j = 1, . . . , M. Then, we 
have 

£ 

y e (s,u/) = ^2{y{s,u.j)^ i )ip l . 
i=l 

We write j/(s, Uj) to stress the dependence on the constant control Uj. Each 
trajectory ?/(s, Uj) has some coefficients (t) for i = 1, . . . ,£, j = 1, . . . , M. 
The coefficients w^(s) will belong to intervals of the type [wj , wp-*] where 
we chose for i = 1, . . . ,£, cti, bi such that: 

- . r (1) (itf)-i 

a,i = mm{w^ ,...,w- j 

Oj = max{w- . . . , '}. 

Then, we have a method to compute the intervals and we turn our attention to 
the numerical solution of an optimal control problem for evolutive equation, 
as we will see in the following section. 



4. Adapting POD approximation 

We now present an adaptive method to compute POD basis. Since our final 
goal is to obtain the optimal feedback law by means of HJB equations, we 
will have a big constraint on the number of variables in the state space for 
numerical solution of an HJB. 

We will see that, for a parabolic equation, one can try to solve the problem 
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with only three/four POD basis functions; they are enough to describe the 
solution in a rather accurate way. In fact the singular values decay pretty 
soon and it's easier to work with a really low-rank dimensional problem. 
On the contrary, hyperbolic equations do not have this nice property for 
their singular values and they will require a rather large set of POD basis 
functions to get accurate results. Note that we can not follow the approach 
suggested in [13] because we can not add more basis functions when it turns 
to be necessary due to the constraint already mentioned. Then, it is quite 
natural to split the problem into subproblems having different POD basis 
functions. The crucial point is to decide the splitting in order to have the same 
number of basis functions in each subdomain with a guaranteed accuracy in 
the approximation. 

Let us first give an illustrative example for the parabolic case, considering a 
ID advection-diffusion equation: 

y s (x, s) - ey xx (x, s) + cy x (x, s) = ^ 
y(x,0) = y (x), 

with x e [a,b],s £ [0,T],e,c G E. 

We use a finite difference approximation for this equation based on an explicit 
Euler method in time combined with the standard centered approximation 
of the second order term and with an up-wind correction for the advection 
term. The snapshots will be taken from the sequence generated by the finite 
difference method. The final time is T — 5, moreover a = — 1, b = 4. The 
initial condition is yo(x) — 5x — 5x 2 , when < x < 1, otherwise. 
For e — 0.05 and c = 1 with only 3 POD basis functions, the approximation 
fails (see Figure [T]) . Note that in this case the advection is dominating the 
diffusion, a low number of POD basis functions will not suffice to get an 
accurate approximation (Figure l.b). However, the adaptive method which 
only uses 3 POD basis functions will give accurate results (Figure l.d). 

The idea which is behind the adaptive method is the following: we do not 
consider all the snapshots together in the whole interval [0, T] but we group 
them. Instead of taking into account the whole interval [0,T], we prefer to 
split it in sub-intervals 

[0,T] = uf =0 [T fc ,T fc+1 ] 

where K is a-priori unknown, Tq — 0, Tk = T and = ti for some i. In this 
way, choosing properly the length of the k— th interval [Tk, Tfc+i], we consider 
only the snapshots falling in that sub-interval, typically there will be at least 
three snapshots in every sub-interval. Then we have enough informations 
in every sub-interval and we can apply the standard routines (explained in 
Section 2) to get a "local" POD basis. 

Now let us explain how to divide our time interval [0,T]. We will choose 
a parameter to check the accuracy of the POD approximation and define a 
threshold. Above that threshold we loose in accuracy and we need to compute 



a new POD basis. A good parameter to check the accuracy is £{t) (see < 2.4 1), 
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(a) 



(b) 





(c) 



(d) 



FIGURE 1. Equation (4.1|:(a) solved with finite difference; 
(b) POD-Galerkin approximation with 3 POD basi; (c) 
solved via POD-Galerkin approximation with 5 POD basis; 
(d) Adapting 3 POD basis functions. 



as it was suggested by several authors. The method to define the splitting of 
[0, T] and the size of every sub-interval works as follows. We start computing 
the SVD of the matrix Y that gives us informations about our dynamics in 
the whole time interval. We check the accuracy at every tj, i = 1, . . . N, and 
if at tk the indicator is above the tolerance we set T\ = tk and we divide the 
interval in two parts, [0, T{) and (T\, T\. Now we just consider the snapshots 
related the solution up to the time T\. Then we iterate this idea until the 
indicator is below the threshold. When the first interval is found, we restart 
the procedure in the interval [T\ , T] and we stop when we reach the final 
time T. Note that the extrema of every interval coincide by construction 
with one of our discrete times t, = iAt so that the global solution is easily 
obtained linking all the sub-problems which always have a snapshot as initial 
condition. A low value for the threshold will also guarantee that we will not 
have big jumps passing from one sub-interval to the next. 
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This idea can be applied also when we have a controlled dynamic (see (5.1 1). 
First of all we have to decide how to collect the snapshots, since the control 
u(t) is completely unknown. One can make a guess and use the dynamics 
and the functional corresponding to that guess, by these informations we 
can compute the POD basis. Once the POD basis is obtained we will get 
the optimal feedback law after having solved a reduced HJB equation as we 
already explained. Let us summarize the method in the following step-by-step 
presentation. 



ALGORITHM 

Start: Inizialization 

Step 1: collect the snapshots in [0,T] 

Step 2: divide [0,T] according to £{£) 

For i=0 to N-l 

Do 

Step 3: apply SVD to get the POD basis in each sub-interval [ti,ij 
Step 4: discretize the space of controls 

Step 5: project the dynamics onto the (reduced) POD space 
Step 6: select the intervals for the POD reduced variables 
Step 7: solve the corresponding HJB in the reduced space 

for the interval [ti,t i+ i] 
Step 8: go back to the original coordinate space 
End 



5. Numerical experiments 

In this section we present some numerical tests for the controlled heat equa- 
tion and for the advection-diffusion equation with a quadratic cost functional. 
Consider the following advection-diffusion equation: 

/ Vs(x,s) - ey xx (x, s) + cy x (x,s) = u(s) , . 

\v(x,0)=Vo{x), [bA) 

with x G [a, b], s G [0, T], e G R+ and ceR. 

Note that changing the parameters c and e we can obtain the heat equation 
(c = 0) and the advection equation (e = 0). The functional to be minimized 
is 

Jy ,tH-))= [ \\y(x,s)~y(x,s)\\ 2 +R\\u( s )\\ 2 ds, (5.2) 
Jo 

i.e. we want to stay close to a reference trajectory y while minimizing the norm 
of u. Note that we dropped the discount factor setting A = 0. Typically in our 
test problems y is obtained by applying a particular control u to the dynamics. 
The numerical simulations reported in this papers have been made on a server 
SUPERMICRO 8045C-3RB with 2 cpu Intel Xeon Quad-Core 2.4 Ghz and 32 
GB RAM under SLURM (https://computing.llnl.gov/linux/slurm/). 



12 



A. Alia and M. Falcone 



Test 1: Heat equation with smooth initial data. We compute the snapshots 
with a centered/forward Euler scheme with space step Ax — 0.02, and time 
step At = 0.012, e = 1/60, c = 0,R= 0.01 and T = 5. The initial condition 
is ya(x) — 5x — 5x 2 , and y{x,s) = 0. In Figure [2] we compare four different 
approximations concerning the heat equation: (a) is the solution for u(t) — 
0, (b) is its approximation via POD (non adaptive), (c) is the direct LQR 
solution computed by MATLAB without POD and, finally, the approximate 
optimal solution obtained coupling POD and HJB. The approximate value 
function is computed for At = 0.1 Ax = 0.1 whereas the optimal trajectory 
as been obtained with At — 0.01. Test 1, and even Test 2, have been solved 
in about half an hour of CPU time. 

Note that in this example the approximate solution is rather accurate because 
the regularity of the solution is high due to the diffusion term. Since in the 
limit the solution tends to the average value the choice of the snapshots 
will not affect too much the solution, i.e. even with a rough choice of the 
snapshots will give us a good approximation. The difference between Figure 
2c and Figure 2d is due to the fact that the control space is continuous for 
2c and discrete for 2d. 

Test 2: Heat equation with no-smooth intial data. In this section we change 
the initial condition with a function which is only Lipschitz continuos: yo{x) 
I — According to Test 1, we consider the same parameters, (see Figure 
§. 

Riccati's equation has been solved by a MATLAB LQR routine. Thus, we 
have used the solution given by this routine as the correct solution in order to 
compare the errors in L 1 and L 2 norm between the reduced Riccati's equa- 
tion and our approach based on the reduced HJB equation. Since we do not 
have any information, the snapshots are computed for u = 0. This is only a 
guess, but in the parabolic case fits well due to the diffusion term. 
As in Test 1, the choice of the snapshots docs not effect strongly the ap- 





L 1 


L 2 


yLQR _ yPOD+LQR 


0.0221 


0.0172 


yLQR _ POD+HJB 


0.0204 


0.0171 



Table 1. Test 2: L 1 and L 2 errors at time T for the optimal 
approximate solution. 



proximation due to the asymptotic behavior of the solution. The presence of 
a Lipschitz continuous initial condition has almost no influence on the global 
error (see Table 1). 
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(a) 



(b) 





(c) 



(d) 



Figure 2. Test l:(a) Heat Equation without control; (b) 
Heat Equation without control, 3 POD basis; (c) Con- 
trolled solution with LQR-MATLAB; (d) Approximate so- 
lution POD (3 basis functions) + HJB. 



Test 3: Advection-Diffusion equation. The advection-diffusion equation needs 
a different method. We can not use the same y we had in the parabolic case, 
mainly because in Riccati's equation the control is free and is not bounded, on 
the contrary when we solve an HJB we have to discretize the space of controls. 
We modified the problem in order to deal with bang-bang controls. We get 



y in (5.2) just plugging in the control u = 0. We have considered the control 
space corresponding only to three values in [—1, 1], then U = {—1, 0, 1}. We 
first have tried to get a controlled solution, without any adaptive method 
and, as expected, we obtained a bad approximation (see Figure [4]). From 
Figure [4] it's clear that POD with four basis functions is not able to catch 
the behavior of the dynamics, so we have applied our adaptive method. 
We have consider: T = 3, Ax = 0.1, At = 0.008, a = -1, b = 4, R = 
0.01. According to our algorithm, the time interval [0,3] was divided into 
[0, 0.744] U [0.744, 1.496] U [1.496, 3]. As we can see our last interval is bigger 
than the others, this is due to the diffusion term (see Figure[5|. The L 2 — error 
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(c) (d) 

Figure 3. Test 2: (a) exact solution for u = 0; (b) Exact 
solution for u — POD (3 basis functions); (c) Approxi- 
mate optimal solution for LQR-MATLAB; (d) Approximate 
solution POD (3 basis functions)+ HJB. 




FIGURE 4. Test 3: Solution y on the left, approximate solu- 
tion on the right with POD (4 basis functions) 
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Figure 5. Test 3: Solution for u 
optimal solution (right). 



(left), approximate 



is 0.0761, and the computation of the optimal solution via HJB has required 
about six hours of CPU time. In Figure 4 we compare the exact solution 
with the numerical solution based on a POD representation. Note that, in 
this case, the choice of only 4 basis functions for the whole interval [0, T] 
gives a very poor result due to the presence of the advection term. Looking 
at Figure 5 one can see the improvement of our adaptive technique which 
takes always 4 basis functions in each sub-interval. 

In order to check the quality of our approximation we have computed the 
numerical residual, defined as: 

= hs{x,s) - ey xx {x, s) + cy x (x,s) - u(s)\\. 

The residual for the solution of the control problem computed without our 
adaptive technique is 1.1, whereas the residual for the adaptive method is 
2 * 10~ 2 . As expected from the pictures, there is a big difference between 
these two value. 

Test 4: Advection-DifFusion equation. In this test we take a different y, 



namely the solution of (5.1 ) corresponding to the control 



-1 < t < 1 
u(t) = { 1 < t < 2 
1 2 < t < 3. 

We want to emphasize we can obtain nice results when the space of controls 
has few element. The parameters were the same used in Test 3. The L 2 — error 
is 0.09, and the time was the same we had in Test 3. In Figure [6] we can see 
our approximation. In Figure 6 one can see that the adaptive technique can 
also deal with discontinuous controls. 



In this test, the residual for the solution of the control problem without 
our adaptive technique is 2, whereas the residual for the adaptive method 
is 3 * 10~ 2 . Again, the residual shows the higher accuracy of the adaptive 
routine. 



1G 
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Figure 6. Test 4: Solution for u (left), approximate optimal 
solution (right). 

6. Conclusions 

As we have discussed, a reasonable coupling between POD and HJB equation 
can produce feedback controls for infinite dimensional problem. For advec- 
tion dominated equations that simple idea has to be implemented in a clever 
way to be successful. It particular, the application of an adaptive technique 
is crucial to obtain accurate approximations with a low number of POD basis 
functions. This is still an essential requirement when dealing with the Dy- 
namic Programming approach, which suffers from the curse-of-dimensionality 
although recent developments in the methods used for HJB equations will al- 
low to increase this bound in the next future (for example by applying patchy 
techniques) . 

Another important point is the discretization of the control space. In 
our examples, the number of optimal control is rather limited and this will be 
enough for problems which have a bang-bang structure for optimal controls. 
In general, we will need also an approximation of the control space via reduced 
basis methods. This point as well as a more detailed analysis of the procedure 
outlined in this paper will be addressed in our future work. 
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